Wolbachia (oligotyping)

Load packages, paths, functions

# Load main packages, paths and custom functions
source("../../../source/main_packages.R")
source("../../../source/paths.R")
source("../../../source/functions.R")

# Load supplementary packages
packages <- c("RColorBrewer", "ggpubr", "cowplot", "Biostrings", "openxlsx", "kableExtra")
invisible(lapply(packages, require, character.only = TRUE))

Preparation

Tables preparation

Seqtab

# move to oligotyping directory
setwd(paste0(path_oligo,"/wolbachia/oligotyping_Wolbachia_sequences-c2-s1-a0.0-A0-M10"))

# load the matrix count table
matrix_count <- read.table("MATRIX-COUNT.txt", header = TRUE) %>% t()

# arrange it
colnames(matrix_count) <- matrix_count[1,]
matrix_count <- matrix_count[-1,]
matrix_count <- matrix_count %>% as.data.frame()

# print it
matrix_count %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
CTC10 CTC11 CTC12 CTC13 CTC14 CTC15 CTC2 CTC3 CTC4 CTC5 CTC6 CTC7 CTC9 NP14 NP2 NP20 NP27 NP29 NP30 NP34 NP35 NP37 NP38 NP39 NP41 NP42 NP43 NP44 NP5 NP8 S100 S101 S102 S104 S105 S106 S107 S108 S109 S110 S113 S121 S122 S123 S124 S126 S127 S128 S146 S147 S148 S150 S151 S152 S153 S154 S160 S162 S163 S164 S165 S166 S167 S169 S170 S18 S19 S20 S21 S22 S23 S24 S25 S26 S27 S28 S29 S30 S31 S32 S33 S34 S35 S36 S37 S38 S39 S40 S41 S42 S43 S44 S45 S46 S47 S48 S49 S50 S51 S52 S53 S55 S56 S57 S58 S59 S60 S61 S63 S64 S65 S79 S80 S83 S84 S85 S86 S87 S88 S89
AT 75 88 16 37 149 164 211 23 9 24 29 2 23 7510 618873 60 1132 9 64 29 637 101 5924 479 3 80 1571 2 703792 312332 52168 4204 3446 52042 55188 32800 51828 55473 58672 57141 11538 16105 7915 16231 14357 12030 19909 13241 20256 20238 11212 12053 18527 12233 13371 7439 57333 20325 9679 18121 14055 10723 10940 9242 7090 4234 30679 39952 22351 9576 35749 38103 47217 5222 5321 4331 18 4200 4782 2624 12821 5495 5240 4692 764 45671 11815 19172 7 4091 9243 7837 17710 498 1850 56207 33144 28127 60995 21063 4 50226 42384 48894 30175 32592 44266 49111 53121 47253 11 59372 52605 41938 56299 41407 61573 65420 52538 4
GT 10 15 4 10 25 41 42 5 2 1 7 0 3 195 15502 0 37 1 2 0 6 1 147 16 0 2 57 0 18486 14077 100 6 2 140 113 88 105 75 129 89 32 3040 1385 2732 2582 2206 3493 2227 3515 3517 2058 2253 3254 2047 2700 1570 10526 3468 1829 3061 2613 2266 2164 1511 1157 6 56 90 35 25 66 87 106 6 10 4 1 2 7 2 2 3 4 8 3 67 21 27 0 9 7 5 36 2 1 39 28 21 129 36 0 88 110 112 71 71 72 117 142 176 0 121 58 129 136 98 130 172 150 0
AG 2 1 1 1 4 9 6 0 5 2 2 0 2 186 11074 4 36 0 0 1 10 2 88 11 1 3 53 0 11433 6724 85 5 3 51 68 64 97 71 95 54 17 502 173 408 435 372 472 285 438 491 314 297 402 245 369 255 1623 592 314 431 454 450 367 352 254 6 35 31 31 12 66 45 55 8 8 4 0 8 10 3 12 7 8 2 2 38 17 22 0 1 5 5 7 0 3 60 39 23 67 29 0 33 25 53 29 46 37 42 32 70 0 49 42 56 63 41 41 69 67 0
TT 4 3 1 0 6 3 6 0 3 0 1 0 1 1 318 0 2 0 0 0 0 0 4 0 0 0 0 0 312 176 6 0 0 1 7 6 8 10 11 6 2 582 288 626 644 527 642 480 678 794 492 409 653 478 522 284 2539 660 418 620 476 420 448 369 292 0 3 7 2 2 7 4 7 0 2 0 0 0 2 0 2 1 1 0 0 4 2 2 0 1 4 1 2 0 0 3 1 4 8 2 0 6 2 7 2 5 6 4 13 9 0 10 3 7 12 8 3 8 4 0
CC 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1080 0 0 0 0 0 3 0 5 1 0 0 7 0 1082 608 73 12 3 116 108 68 87 84 93 72 21 18 9 47 42 26 48 18 30 31 18 16 22 24 18 18 130 31 15 38 31 22 19 12 16 5 63 84 31 17 51 99 109 8 9 4 0 3 7 3 19 4 3 3 3 62 21 32 0 4 20 12 41 0 2 69 35 16 125 52 0 75 81 74 58 57 78 236 97 83 0 138 74 87 113 76 152 177 170 0
GG 1 0 0 0 0 1 2 0 0 1 0 0 0 9 681 0 1 0 0 0 3 0 6 0 0 0 3 0 804 646 0 0 0 0 1 0 0 0 0 0 0 99 31 74 72 68 114 47 71 80 55 49 78 47 55 49 310 91 76 76 85 82 97 65 43 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 0 0 0 1 0 0
AC 0 0 0 0 0 0 0 0 0 0 1 0 0 0 411 0 0 0 0 0 0 0 2 0 0 0 3 0 248 216 54 0 2 49 92 27 27 22 22 14 5 15 2 11 12 2 18 7 8 18 6 4 7 5 3 5 37 11 1 8 8 8 8 1 0 0 19 28 6 10 24 35 62 6 5 1 0 1 2 0 0 0 1 5 0 23 1 18 0 0 0 2 21 0 0 13 2 2 73 13 0 19 7 14 21 27 25 35 39 37 0 64 5 55 53 60 84 111 75 0

Taxonomy

# move to oligotyping directory
setwd(paste0(path_oligo,"/wolbachia/oligotyping_Wolbachia_sequences-c2-s1-a0.0-A0-M10"))

# load the fasta table
fasta <- readDNAStringSet("OLIGO-REPRESENTATIVES.fasta")

# arrange it
fasta <- fasta %>% as.data.frame()
colnames(fasta) <- "seq"
fasta$oligotype <- rownames(fasta)
fasta <- fasta %>% dplyr::select(-c(seq))

# print it
fasta %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
oligotype
AT AT
GT GT
AG AG
TT TT
CC CC
GG GG
AC AC

Change oligotype name by oligotype / MED nodes in the matrix count

# Reference file 

## move to tsv directory
setwd(path_tsv)

## load the reference table
ref_oligo_med2 <- read.table("2B_REF_info_wolbachia.tsv", sep="\t", header = TRUE)

## select only the 7 oligotypes of Wolbachia
ref_oligo_med2 <- ref_oligo_med2[!is.na(ref_oligo_med2$oligotype),]

## change order of columns
ref_oligo_med2 <- ref_oligo_med2 %>% select(c(seq, oligotype, MED_node_frequency_size, OLIGO_oligotype_frequency_size))

## create a column with reference name (will be used in plots)
ref_oligo_med2$ref <- paste0("oligotype_", ref_oligo_med2$OLIGO_oligotype_frequency_size, " / node_", ref_oligo_med2$MED_node_frequency_size)

## create a copy of fasta 
fasta2 <- fasta

# Matrix count

## create an oligotype column in the matrix count
matrix_count$oligotype <- rownames(matrix_count)

## change order of columns
matrix_count <- matrix_count %>% dplyr::select(c(oligotype, everything()))

## merge the matrix count and the reference dataframe
matrix_count2 <- matrix_count %>% merge(ref_oligo_med2 %>% dplyr::select(-c(seq)), by="oligotype")

## change order of columns
matrix_count2 <- matrix_count2 %>% dplyr::select(c(oligotype, MED_node_frequency_size, OLIGO_oligotype_frequency_size, ref, everything()))

## change rownames
rownames(matrix_count2) <- matrix_count2$ref

## change order of columns
matrix_count2 <- matrix_count2 %>% dplyr::select(-c(oligotype, ref, MED_node_frequency_size, OLIGO_oligotype_frequency_size))

## print it
matrix_count2 %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
CTC10 CTC11 CTC12 CTC13 CTC14 CTC15 CTC2 CTC3 CTC4 CTC5 CTC6 CTC7 CTC9 NP14 NP2 NP20 NP27 NP29 NP30 NP34 NP35 NP37 NP38 NP39 NP41 NP42 NP43 NP44 NP5 NP8 S100 S101 S102 S104 S105 S106 S107 S108 S109 S110 S113 S121 S122 S123 S124 S126 S127 S128 S146 S147 S148 S150 S151 S152 S153 S154 S160 S162 S163 S164 S165 S166 S167 S169 S170 S18 S19 S20 S21 S22 S23 S24 S25 S26 S27 S28 S29 S30 S31 S32 S33 S34 S35 S36 S37 S38 S39 S40 S41 S42 S43 S44 S45 S46 S47 S48 S49 S50 S51 S52 S53 S55 S56 S57 S58 S59 S60 S61 S63 S64 S65 S79 S80 S83 S84 S85 S86 S87 S88 S89
oligotype_AC (80) | size:2504 / node_N0253 (80) | size:2481 0 0 0 0 0 0 0 0 0 0 1 0 0 0 411 0 0 0 0 0 0 0 2 0 0 0 3 0 248 216 54 0 2 49 92 27 27 22 22 14 5 15 2 11 12 2 18 7 8 18 6 4 7 5 3 5 37 11 1 8 8 8 8 1 0 0 19 28 6 10 24 35 62 6 5 1 0 1 2 0 0 0 1 5 0 23 1 18 0 0 0 2 21 0 0 13 2 2 73 13 0 19 7 14 21 27 25 35 39 37 0 64 5 55 53 60 84 111 75 0
oligotype_AG (109) | size:42030 / node_N0245 (109) | size:41133 2 1 1 1 4 9 6 0 5 2 2 0 2 186 11074 4 36 0 0 1 10 2 88 11 1 3 53 0 11433 6724 85 5 3 51 68 64 97 71 95 54 17 502 173 408 435 372 472 285 438 491 314 297 402 245 369 255 1623 592 314 431 454 450 367 352 254 6 35 31 31 12 66 45 55 8 8 4 0 8 10 3 12 7 8 2 2 38 17 22 0 1 5 5 7 0 3 60 39 23 67 29 0 33 25 53 29 46 37 42 32 70 0 49 42 56 63 41 41 69 67 0
oligotype_AT (120) | size:3890567 / node_N0711 (120) | size:3744518 75 88 16 37 149 164 211 23 9 24 29 2 23 7510 618873 60 1132 9 64 29 637 101 5924 479 3 80 1571 2 703792 312332 52168 4204 3446 52042 55188 32800 51828 55473 58672 57141 11538 16105 7915 16231 14357 12030 19909 13241 20256 20238 11212 12053 18527 12233 13371 7439 57333 20325 9679 18121 14055 10723 10940 9242 7090 4234 30679 39952 22351 9576 35749 38103 47217 5222 5321 4331 18 4200 4782 2624 12821 5495 5240 4692 764 45671 11815 19172 7 4091 9243 7837 17710 498 1850 56207 33144 28127 60995 21063 4 50226 42384 48894 30175 32592 44266 49111 53121 47253 11 59372 52605 41938 56299 41407 61573 65420 52538 4
oligotype_CC (92) | size:7065 / node_N0026 (92) | size:7069 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1080 0 0 0 0 0 3 0 5 1 0 0 7 0 1082 608 73 12 3 116 108 68 87 84 93 72 21 18 9 47 42 26 48 18 30 31 18 16 22 24 18 18 130 31 15 38 31 22 19 12 16 5 63 84 31 17 51 99 109 8 9 4 0 3 7 3 19 4 3 3 3 62 21 32 0 4 20 12 41 0 2 69 35 16 125 52 0 75 81 74 58 57 78 236 97 83 0 138 74 87 113 76 152 177 170 0
oligotype_GG (44) | size:4080 / node_N0250 (44) | size:4011 1 0 0 0 0 1 2 0 0 1 0 0 0 9 681 0 1 0 0 0 3 0 6 0 0 0 3 0 804 646 0 0 0 0 1 0 0 0 0 0 0 99 31 74 72 68 114 47 71 80 55 49 78 47 55 49 310 91 76 76 85 82 97 65 43 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 0 0 0 1 0 0
oligotype_GT (111) | size:119651 / node_N0241 (111) | size:114576 10 15 4 10 25 41 42 5 2 1 7 0 3 195 15502 0 37 1 2 0 6 1 147 16 0 2 57 0 18486 14077 100 6 2 140 113 88 105 75 129 89 32 3040 1385 2732 2582 2206 3493 2227 3515 3517 2058 2253 3254 2047 2700 1570 10526 3468 1829 3061 2613 2266 2164 1511 1157 6 56 90 35 25 66 87 106 6 10 4 1 2 7 2 2 3 4 8 3 67 21 27 0 9 7 5 36 2 1 39 28 21 129 36 0 88 110 112 71 71 72 117 142 176 0 121 58 129 136 98 130 172 150 0
oligotype_TT (89) | size:15422 / node_N0244 (89) | size:15236 4 3 1 0 6 3 6 0 3 0 1 0 1 1 318 0 2 0 0 0 0 0 4 0 0 0 0 0 312 176 6 0 0 1 7 6 8 10 11 6 2 582 288 626 644 527 642 480 678 794 492 409 653 478 522 284 2539 660 418 620 476 420 448 369 292 0 3 7 2 2 7 4 7 0 2 0 0 0 2 0 2 1 1 0 0 4 2 2 0 1 4 1 2 0 0 3 1 4 8 2 0 6 2 7 2 5 6 4 13 9 0 10 3 7 12 8 3 8 4 0 
## edit the fasta dataframe
fasta2 <- fasta2 %>% merge(ref_oligo_med2 %>% dplyr::select(-c(seq)), by="oligotype")
rownames(fasta2) <- fasta2$ref
fasta2 <- fasta2 %>% dplyr::select(-c(MED_node_frequency_size, OLIGO_oligotype_frequency_size, oligotype))

## print it
fasta2 %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
ref
oligotype_AC (80) | size:2504 / node_N0253 (80) | size:2481 oligotype_AC (80) | size:2504 / node_N0253 (80) | size:2481
oligotype_AG (109) | size:42030 / node_N0245 (109) | size:41133 oligotype_AG (109) | size:42030 / node_N0245 (109) | size:41133
oligotype_AT (120) | size:3890567 / node_N0711 (120) | size:3744518 oligotype_AT (120) | size:3890567 / node_N0711 (120) | size:3744518
oligotype_CC (92) | size:7065 / node_N0026 (92) | size:7069 oligotype_CC (92) | size:7065 / node_N0026 (92) | size:7069
oligotype_GG (44) | size:4080 / node_N0250 (44) | size:4011 oligotype_GG (44) | size:4080 / node_N0250 (44) | size:4011
oligotype_GT (111) | size:119651 / node_N0241 (111) | size:114576 oligotype_GT (111) | size:119651 / node_N0241 (111) | size:114576
oligotype_TT (89) | size:15422 / node_N0244 (89) | size:15236 oligotype_TT (89) | size:15422 / node_N0244 (89) | size:15236

Metadata

metadata <- read.csv(paste0(path_metadata,"/metadata_08_02_2021.csv"), sep=";")
rownames(metadata) <- metadata$Sample

Phyloseq object with oligotypes

# convert matrix_count into matrix and numeric
matrix_count <- matrix_count2 %>% as.matrix()
class(matrix_count) <- "numeric"

# phyloseq elements
OTU = otu_table(as.matrix(matrix_count), taxa_are_rows =TRUE)
TAX = tax_table(as.matrix(fasta2))
SAM = sample_data(metadata)

# phyloseq object
ps <- phyloseq(OTU, TAX, SAM)
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 7 taxa and 120 samples ]
## sample_data() Sample Data:       [ 120 samples by 15 sample variables ]
## tax_table()   Taxonomy Table:    [ 7 taxa by 1 taxonomic ranks ]
compute_read_counts(ps)
## [1] 4081319
# remove blanks
ps <- subset_samples(ps, Location!="Blank")
ps <- check_ps(ps)
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 7 taxa and 113 samples ]
## sample_data() Sample Data:       [ 113 samples by 15 sample variables ]
## tax_table()   Taxonomy Table:    [ 7 taxa by 1 taxonomic ranks ]

Create new metadata with Percent

Load ps with all samples (for final plot)

setwd(path_rdata)
ps.filter <- readRDS("1D_MED_phyloseq_decontam.rds")
ps.filter <- check_ps(ps.filter)

Edit new metadata with Percent_wolbachia

guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 16, face = "italic", colour = "Black", angle = 0)))

# add read depth in sample table of phyloseq object
sample_data(ps.filter)$Read_depth <- sample_sums(ps.filter)

# select Wolbachia
ps.wolbachia <- ps.filter %>% subset_taxa(Genus=="Wolbachia")

# add read depth of Wolbachia
sample_data(ps.filter)$Read_wolbachia <- sample_sums(ps.wolbachia)
sample_data(ps.filter) %>% colnames()
##  [1] "Sample"         "Well"           "Primer1"        "Primer2"       
##  [5] "Location"       "Field"          "Country"        "Organ"         
##  [9] "Species"        "Individual"     "Individuals"    "Date"          
## [13] "Run"            "Control"        "Dna"            "Read_depth"    
## [17] "Read_wolbachia"
sample_data(ps.wolbachia) %>% colnames()
##  [1] "Sample"      "Well"        "Primer1"     "Primer2"     "Location"   
##  [6] "Field"       "Country"     "Organ"       "Species"     "Individual" 
## [11] "Individuals" "Date"        "Run"         "Control"     "Dna"        
## [16] "Read_depth"
# add percent of Wolbachia
sample_data(ps.filter)$Percent_wolbachia <- sample_data(ps.filter)$Read_wolbachia / sample_data(ps.filter)$Read_depth

# round the percent of Wolbachia at 2 decimals
sample_data(ps.filter)$Percent_wolbachia <- sample_data(ps.filter)$Percent_wolbachia %>% round(2)

# extract metadata table
test <- data.frame(sample_data(ps.filter))

# merge this metadata table with the other
new.metadata <- data.frame(sample_data(ps)) %>% merge(test %>% dplyr::select(c(Sample, Read_depth, Read_wolbachia, Percent_wolbachia)), by="Sample")
new.metadata <- test[new.metadata$Sample %in% sample_names(ps),]
rownames(new.metadata) <- new.metadata$Sample

# print it
new.metadata %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
Sample Well Primer1 Primer2 Location Field Country Organ Species Individual Individuals Date Run Control Dna Read_depth Read_wolbachia Percent_wolbachia
CTC1 CTC1 G5 V4-SA707 V3-SA505 Wolbachia - Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 58 29779 0 0.00
CTC10 CTC10 D6 V4-SA704 V3-SA506 Wolbachia - Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 37,9 2609 92 0.04
CTC11 CTC11 E6 V4-SA705 V3-SA506 Wolbachia - Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 58 13874 107 0.01
CTC12 CTC12 F6 V4-SA706 V3-SA506 Wolbachia - Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 40,1 1146 22 0.02
CTC13 CTC13 G6 V4-SA707 V3-SA506 Wolbachia - Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 58 18035 48 0.00
CTC14 CTC14 H6 V4-SA708 V3-SA506 Wolbachia - Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 6,17 1708 184 0.11
CTC15 CTC15 I6 V4-SA709 V3-SA506 Wolbachia - Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 58 23180 218 0.01
CTC2 CTC2 H5 V4-SA708 V3-SA505 Wolbachia - Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 58 30692 268 0.01
CTC3 CTC3 I5 V4-SA709 V3-SA505 Wolbachia - Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 58 39920 28 0.00
CTC4 CTC4 J5 V4-SA710 V3-SA505 Wolbachia - Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 1,15 2139 19 0.01
CTC5 CTC5 K5 V4-SA711 V3-SA505 Wolbachia - Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 58 15789 28 0.00
CTC6 CTC6 L5 V4-SA712 V3-SA505 Wolbachia - Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 58 19753 40 0.00
CTC9 CTC9 C6 V4-SA703 V3-SA506 Wolbachia - Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 33,6 4980 29 0.01
NP14 NP14 K4 V4-SA711 V3-SA504 Guadeloupe Field Guadeloupe Ovary Aedes aegypti 1a 0 N/A run3 True sample 0,437 7973 7901 0.99
NP2 NP2 K3 V4-SA711 V3-SA503 Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus 1c 0 N/A run3 True sample 41,6 648335 647941 1.00
NP20 NP20 E5 V4-SA705 V3-SA505 Guadeloupe Field Guadeloupe Ovary Aedes aegypti 3a 0 N/A run3 True sample 0,357 136 64 0.47
NP27 NP27 L5 V4-SA712 V3-SA505 Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus 7c 0 N/A run3 True sample 1,16 1234 1208 0.98
NP29 NP29 B6 V4-SA702 V3-SA506 Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus 9c 0 N/A run3 True sample 0,314 203 10 0.05
NP30 NP30 C6 V4-SA703 V3-SA506 Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus 10c 0 N/A run3 True sample 0,666 228 66 0.29
NP34 NP34 G6 V4-SA707 V3-SA506 Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus 14c 0 N/A run3 True sample 0,486 95 30 0.32
NP35 NP35 H6 V4-SA708 V3-SA506 Guadeloupe Field Guadeloupe Whole Aedes aegypti 7a 0 N/A run3 True sample 4,64 196532 659 0.00
NP36 NP36 I6 V4-SA709 V3-SA506 Guadeloupe Field Guadeloupe Whole Aedes aegypti 8a 0 N/A run3 True sample 1,06 249 0 0.00
NP37 NP37 J6 V4-SA710 V3-SA506 Guadeloupe Field Guadeloupe Whole Aedes aegypti 9a 0 N/A run3 True sample 22,7 419340 104 0.00
NP38 NP38 K6 V4-SA711 V3-SA506 Guadeloupe Field Guadeloupe Whole Aedes aegypti 10a 0 N/A run3 True sample 3,88 282479 6176 0.02
NP39 NP39 L6 V4-SA712 V3-SA506 Guadeloupe Field Guadeloupe Whole Aedes aegypti 11a 0 N/A run3 True sample 20,2 218684 507 0.00
NP41 NP41 B7 V4-SA702 V3-SA507 Guadeloupe Field Guadeloupe Whole Aedes aegypti 13a 0 N/A run3 True sample 5,32 247152 4 0.00
NP42 NP42 C7 V4-SA703 V3-SA507 Guadeloupe Field Guadeloupe Whole Aedes aegypti 14a 0 N/A run3 True sample 4,65 185157 85 0.00
NP43 NP43 D7 V4-SA704 V3-SA507 Guadeloupe Field Guadeloupe Whole Aedes aegypti 15a 0 N/A run3 True sample 6,89 239335 1694 0.01
NP44 NP44 E7 V4-SA705 V3-SA507 Guadeloupe Field Guadeloupe Whole Aedes aegypti 16a 0 N/A run3 True sample 21,7 156879 2 0.00
NP5 NP5 B4 V4-SA702 V3-SA504 Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus 2c 0 N/A run3 True sample 33,5 736159 736157 1.00
NP8 NP8 E4 V4-SA705 V3-SA504 Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus 3c 0 N/A run3 True sample 46 334799 334782 1.00
S100 S100 K7 V4-SA711 V3-SA507 Camping Europe Field France Ovary Culex pipiens GL1 1 30/05/2017 run1 True sample 8,02 52486 52486 1.00
S102 S102 A8 V4-SA701 V3-SA508 Camping Europe Field France Ovary Culex pipiens GL2 2 30/05/2017 run1 True sample 0,241 3456 3456 1.00
S104 S104 C8 V4-SA703 V3-SA508 Camping Europe Field France Ovary Culex pipiens GL5 5 30/05/2017 run1 True sample 24,1 52403 52399 1.00
S105 S105 D8 V4-SA704 V3-SA508 Camping Europe Field France Ovary Culex pipiens GL6 6 30/05/2017 run1 True sample 6,83 55577 55577 1.00
S106 S106 E8 V4-SA705 V3-SA508 Camping Europe Field France Ovary Culex pipiens GL7 7 30/05/2017 run1 True sample 51 33053 33053 1.00
S107 S107 F8 V4-SA706 V3-SA508 Camping Europe Field France Ovary Culex pipiens GL8 8 30/05/2017 run1 True sample 32,6 52154 52153 1.00
S108 S108 G8 V4-SA707 V3-SA508 Camping Europe Field France Ovary Culex pipiens GL9 9 30/05/2017 run1 True sample 32,2 55735 55735 1.00
S109 S109 H8 V4-SA708 V3-SA508 Camping Europe Field France Ovary Culex pipiens GL10 10 30/05/2017 run1 True sample 26,5 59023 59023 1.00
S110 S110 I8 V4-SA709 V3-SA508 Camping Europe Field France Ovary Culex pipiens GL11 0 30/05/2017 run1 True sample 27,5 57377 57377 1.00
S121 S121 H1 V4-SA708 V3-SA501 Bosc Field France Ovary Culex pipiens J26 22 28/06/2017 run2 True sample 12,7 20361 20361 1.00
S122 S122 I1 V4-SA709 V3-SA501 Bosc Field France Ovary Culex pipiens J27 23 28/06/2017 run2 True sample 22,7 9803 9803 1.00
S123 S123 J1 V4-SA710 V3-SA501 Bosc Field France Ovary Culex pipiens J28 24 28/06/2017 run2 True sample 6,41 20130 20130 1.00
S124 S124 K1 V4-SA711 V3-SA501 Bosc Field France Ovary Culex pipiens J29 25 28/06/2017 run2 True sample 33,9 18146 18145 1.00
S126 S126 K6 V4-SA711 V3-SA506 Bosc Field France Ovary Culex pipiens J30 26 28/06/2017 run2 True sample 58 15235 15231 1.00
S127 S127 B2 V4-SA702 V3-SA502 Bosc Field France Ovary Culex pipiens J31 27 28/06/2017 run2 True sample 12,8 24696 24696 1.00
S128 S128 C2 V4-SA703 V3-SA502 Bosc Field France Ovary Culex pipiens J32 28 28/06/2017 run2 True sample 35,1 16305 16305 1.00
S146 S146 I3 V4-SA709 V3-SA503 Lavar (labo) Lab France Ovary Culex pipiens MW52 29 29/08/2017 run2 True sample 35,2 25012 24996 1.00
S147 S147 J3 V4-SA710 V3-SA503 Lavar (labo) Lab France Ovary Culex pipiens MW53 30 29/08/2017 run2 True sample 27,1 25171 25169 1.00
S148 S148 K3 V4-SA711 V3-SA503 Lavar (labo) Lab France Ovary Culex pipiens MW54 31 29/08/2017 run2 True sample 43,2 14164 14155 1.00
S150 S150 A4 V4-SA701 V3-SA504 Lavar (labo) Lab France Ovary Culex pipiens MW55 32 29/08/2017 run2 True sample 2,3 15081 15081 1.00
S151 S151 B4 V4-SA702 V3-SA504 Lavar (labo) Lab France Ovary Culex pipiens MW56 33 29/08/2017 run2 True sample 38,8 22944 22944 1.00
S152 S152 C4 V4-SA703 V3-SA504 Lavar (labo) Lab France Ovary Culex pipiens MW57 34 29/08/2017 run2 True sample 39,8 15082 15080 1.00
S153 S153 D4 V4-SA704 V3-SA504 Lavar (labo) Lab France Ovary Culex pipiens MW58 35 29/08/2017 run2 True sample 52 17040 17038 1.00
S154 S154 E4 V4-SA705 V3-SA504 Lavar (labo) Lab France Ovary Culex pipiens MW59 36 29/08/2017 run2 True sample 37,7 9626 9620 1.00
S160 S160 K4 V4-SA711 V3-SA504 Lavar (labo) Lab France Ovary Culex pipiens MW60 37 29/08/2017 run2 True sample 58 72508 72501 1.00
S162 S162 B5 V4-SA702 V3-SA505 Lavar (labo) Lab France Ovary Culex pipiens MW61 38 29/08/2017 run2 True sample 42 25180 25179 1.00
S163 S163 L6 V4-SA712 V3-SA506 Lavar (labo) Lab France Ovary Culex pipiens MW62 39 30/08/2017 run2 True sample 51 12333 12332 1.00
S164 S164 C5 V4-SA703 V3-SA505 Lavar (labo) Lab France Ovary Culex pipiens MW63 40 30/08/2017 run2 True sample 36,6 22368 22355 1.00
S165 S165 D5 V4-SA704 V3-SA505 Lavar (labo) Lab France Ovary Culex pipiens MW64 41 30/08/2017 run2 True sample 53 17731 17723 1.00
S166 S166 E5 V4-SA705 V3-SA505 Camping Europe Field France Ovary Culex pipiens GL4 4 30/05/2017 run2 True sample 23,1 13979 13971 1.00
S167 S167 F5 V4-SA706 V3-SA505 Bosc Field France Ovary Culex pipiens J32 28 28/06/2017 run2 True sample 29,1 14048 14045 1.00
S169 S169 B7 V4-SA702 V3-SA507 Camping Europe Field France Ovary Culex pipiens 5 43 16/05/2017 run2 True sample 5,84 11553 11553 1.00
S170 S170 C7 V4-SA703 V3-SA507 Camping Europe Field France Ovary Culex pipiens 6 44 16/05/2017 run2 True sample 5,55 8852 8852 1.00
S18 S18 A1 V4-SA701 V3-SA501 Lavar (labo) Lab France Whole Culex pipiens MW75 0 30/08/2017 run1 True sample 0,089 4290 4251 0.99
S19 S19 B1 V4-SA702 V3-SA501 Lavar (labo) Lab France Whole Culex pipiens MW65 0 30/08/2017 run1 True sample 21,9 44527 30855 0.69
S20 S20 C1 V4-SA703 V3-SA501 Lavar (labo) Lab France Whole Culex pipiens MW66 0 30/08/2017 run1 True sample 16,6 42864 40192 0.94
S21 S21 D1 V4-SA704 V3-SA501 Lavar (labo) Lab France Whole Culex pipiens MW67 0 30/08/2017 run1 True sample 12,4 33798 22456 0.66
S22 S22 E1 V4-SA705 V3-SA501 Lavar (labo) Lab France Whole Culex pipiens MW68 0 30/08/2017 run1 True sample 24,1 19044 9642 0.51
S23 S23 F1 V4-SA706 V3-SA501 Lavar (labo) Lab France Whole Culex pipiens MW69 0 30/08/2017 run1 True sample 20,8 38172 35964 0.94
S24 S24 G1 V4-SA707 V3-SA501 Lavar (labo) Lab France Whole Culex pipiens MW70 0 30/08/2017 run1 True sample 34,2 42355 38373 0.91
S25 S25 H1 V4-SA708 V3-SA501 Lavar (labo) Lab France Whole Culex pipiens MW71 0 30/08/2017 run1 True sample 21,9 47688 47556 1.00
S26 S26 I1 V4-SA709 V3-SA501 Lavar (labo) Lab France Whole Culex pipiens MW72 0 30/08/2017 run1 True sample 0,322 5394 5250 0.97
S27 S27 J1 V4-SA710 V3-SA501 Lavar (labo) Lab France Whole Culex pipiens MW73 0 30/08/2017 run1 True sample 11,3 24558 5355 0.22
S28 S28 A2 V4-SA701 V3-SA502 Lavar (labo) Lab France Whole Culex pipiens MW74 0 30/08/2017 run1 True sample 0,112 4503 4344 0.96
S30 S30 K1 V4-SA711 V3-SA501 Lavar (labo) Lab France Whole Culex pipiens MW1 0 23/08/2017 run1 True sample 4,43 25353 4214 0.17
S31 S31 L1 V4-SA712 V3-SA501 Lavar (labo) Lab France Whole Culex pipiens MW2 0 23/08/2017 run1 True sample 2,66 20417 4810 0.24
S32 S32 C2 V4-SA703 V3-SA502 Lavar (labo) Lab France Whole Culex pipiens MW3 0 23/08/2017 run1 True sample 0,504 12441 2632 0.21
S33 S33 D2 V4-SA704 V3-SA502 Lavar (labo) Lab France Whole Culex pipiens MW4 0 23/08/2017 run1 True sample 0,782 33867 12856 0.38
S34 S34 E2 V4-SA705 V3-SA502 Lavar (labo) Lab France Whole Culex pipiens MW5 0 23/08/2017 run1 True sample 1,38 9367 5510 0.59
S35 S35 F2 V4-SA706 V3-SA502 Lavar (labo) Lab France Whole Culex pipiens MW6 0 23/08/2017 run1 True sample 0,56 11663 5257 0.45
S36 S36 G2 V4-SA707 V3-SA502 Lavar (labo) Lab France Whole Culex pipiens MW7 0 23/08/2017 run1 True sample 39,8 33020 4710 0.14
S37 S37 H2 V4-SA708 V3-SA502 Lavar (labo) Lab France Whole Culex pipiens MW8 0 23/08/2017 run1 True sample 41,3 18340 772 0.04
S38 S38 I2 V4-SA709 V3-SA502 Lavar (labo) Lab France Whole Culex pipiens MW9 0 23/08/2017 run1 True sample 32,1 54790 45865 0.84
S39 S39 J2 V4-SA710 V3-SA502 Lavar (labo) Lab France Whole Culex pipiens MW10 0 23/08/2017 run1 True sample 37,3 36273 11878 0.33
S40 S40 K2 V4-SA711 V3-SA502 Lavar (labo) Lab France Whole Culex pipiens MW11 0 23/08/2017 run1 True sample 4,58 44448 19274 0.43
S42 S42 A3 V4-SA701 V3-SA503 Camping Europe Field France Whole Culex pipiens GLE1 0 30/05/2017 run1 True sample 0 4107 4106 1.00
S43 S43 B3 V4-SA702 V3-SA503 Camping Europe Field France Whole Culex pipiens GLE2 0 30/05/2017 run1 True sample 0,191 9279 9279 1.00
S44 S44 C3 V4-SA703 V3-SA503 Camping Europe Field France Whole Culex pipiens GLE3 0 30/05/2017 run1 True sample 0,102 8026 7862 0.98
S45 S45 D3 V4-SA704 V3-SA503 Camping Europe Field France Whole Culex pipiens GLE4 0 30/05/2017 run1 True sample 0,223 18150 17817 0.98
S47 S47 F3 V4-SA706 V3-SA503 Camping Europe Field France Whole Culex pipiens GLE6 0 30/05/2017 run1 True sample 0,291 1951 1856 0.95
S48 S48 G3 V4-SA707 V3-SA503 Camping Europe Field France Whole Culex pipiens GLE7 0 30/05/2017 run1 True sample 3,44 56738 56391 0.99
S49 S49 H3 V4-SA708 V3-SA503 Bosc Field France Whole Culex pipiens E1 0 28/06/2017 run1 True sample 1,1 33498 33249 0.99
S50 S50 I3 V4-SA709 V3-SA503 Bosc Field France Whole Culex pipiens E2 0 28/06/2017 run1 True sample 0,771 28481 28193 0.99
S51 S51 J3 V4-SA710 V3-SA503 Bosc Field France Whole Culex pipiens E3 0 28/06/2017 run1 True sample 17,8 61788 61398 0.99
S52 S52 K3 V4-SA711 V3-SA503 Bosc Field France Whole Culex pipiens E4 0 28/06/2017 run1 True sample 0,495 21553 21195 0.98
S55 S55 B4 V4-SA702 V3-SA504 Bosc Field France Whole Culex pipiens E6 0 28/06/2017 run1 True sample 2,85 50447 50447 1.00
S56 S56 C4 V4-SA703 V3-SA504 Bosc Field France Whole Culex pipiens E7 0 28/06/2017 run1 True sample 3,6 42609 42609 1.00
S57 S57 D4 V4-SA704 V3-SA504 Bosc Field France Whole Culex pipiens E8 0 28/06/2017 run1 True sample 4,92 49157 49154 1.00
S58 S58 E4 V4-SA705 V3-SA504 Bosc Field France Whole Culex pipiens E9 0 28/06/2017 run1 True sample 1,63 30357 30357 1.00
S59 S59 F4 V4-SA706 V3-SA504 Bosc Field France Whole Culex pipiens E10 0 28/06/2017 run1 True sample 1,64 32798 32798 1.00
S60 S60 G4 V4-SA707 V3-SA504 Bosc Field France Whole Culex pipiens E11 0 28/06/2017 run1 True sample 2,7 44485 44485 1.00
S61 S61 H4 V4-SA708 V3-SA504 Bosc Field France Whole Culex pipiens E12 0 28/06/2017 run1 True sample 2 49545 49545 1.00
S63 S63 J4 V4-SA710 V3-SA504 Bosc Field France Whole Culex pipiens E14 0 28/06/2017 run1 True sample 6,13 53444 53444 1.00
S64 S64 K4 V4-SA711 V3-SA504 Bosc Field France Whole Culex pipiens E15 0 28/06/2017 run1 True sample 4,15 47628 47628 1.00
S79 S79 B6 V4-SA702 V3-SA506 Camping Europe Field France Ovary Culex pipiens J16 12 28/06/2017 run1 True sample 33,8 59755 59755 1.00
S80 S80 C6 V4-SA703 V3-SA506 Camping Europe Field France Ovary Culex pipiens J17 13 28/06/2017 run1 True sample 4,58 52788 52788 1.00
S83 S83 F6 V4-SA706 V3-SA506 Camping Europe Field France Ovary Culex pipiens J20 16 28/06/2017 run1 True sample 35,5 42272 42272 1.00
S84 S84 G6 V4-SA707 V3-SA506 Camping Europe Field France Ovary Culex pipiens J21 17 28/06/2017 run1 True sample 21 56676 56676 1.00
S85 S85 H6 V4-SA708 V3-SA506 Camping Europe Field France Ovary Culex pipiens J22 18 28/06/2017 run1 True sample 11,6 41690 41690 1.00
S86 S86 I6 V4-SA709 V3-SA506 Camping Europe Field France Ovary Culex pipiens J23 19 28/06/2017 run1 True sample 4,14 61984 61984 1.00
S87 S87 J6 V4-SA710 V3-SA506 Bosc Field France Ovary Culex pipiens J24 20 28/06/2017 run1 True sample 28,1 65958 65958 1.00
S88 S88 K6 V4-SA711 V3-SA506 Bosc Field France Ovary Culex pipiens J25 21 28/06/2017 run1 True sample 8,6 53102 53004 1.00 
# replace metadata in the created phyloseq object
sample_data(ps) <- sample_data(new.metadata)

Taxonomic structure

Count

col <- brewer.pal(7, "Pastel2")

# reshape data for plot
test3 <- test %>% select(c(Sample, Species, Location, Organ, Read_depth, Read_wolbachia)) %>% reshape2::melt(id.vars=c("Sample", "Species", "Location", "Organ"), vars=c("Read_depth", "Read_wolbachia"))

count_whole <- test3[test3$Organ=="Whole",]
count_ovary <- test3[test3$Organ=="Ovary",]

make.italic <- function(x) as.expression(lapply(x, function(y) bquote(italic(.(y)))))

levels(count_whole$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(count_ovary$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(count_whole$Location) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))

levels(count_ovary$Location) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))


# plot
p_count1 <- ggplot(count_whole, aes(x = Sample, y = value, fill=variable))+ 
  geom_bar(position = "dodge", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=12, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text=element_text(size=14), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Location+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
  ylim(0, 900000)+
  geom_text(aes(label=value), position=position_dodge(width=1.1), width=0.25, size=4, hjust=-0.25, vjust=0.5, angle=90)+
  guides(fill=guide_legend(title="Read"))
## Warning: Ignoring unknown parameters: width
p_count2 <- ggplot(count_ovary, aes(x = Sample, y = value, fill=variable))+ 
  geom_bar(position = "dodge", stat = "identity")+
    scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text=element_text(size=14), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
 facet_wrap(~Species+Location+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
    ylim(0, 900000)+
  geom_text(aes(label=value), position=position_dodge(width=0.8), width=0.25, size=4, hjust=-0.25, vjust=0.5, angle=90)+
  guides(fill=guide_legend(title="Read"))
## Warning: Ignoring unknown parameters: width
# afficher plot
p_count1
## Warning: position_dodge requires non-overlapping x intervals
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

p_count2

# panels
p_group <- plot_grid(p_count1+theme(legend.position="none"), 
          p_count2+theme(legend.position="none"), 
          nrow=2, 
          ncol=1)+
    draw_plot_label(c("B1", "B2"), c(0, 0), c(1, 0.5), size = 20)
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals
legend_plot <- get_legend(p_count1 + theme(legend.position="bottom"))
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals
p_counts <- plot_grid(p_group, legend_plot, nrow=2, ncol=1, rel_heights = c(1, .1))
p_counts

Whole (the most abundant nodes)

#fig.height = 12, fig.width=18
guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 16, face = "italic", colour = "Black", angle = 0),
                                            nrow=4, byrow=TRUE))

# select whole
ps.filter.whole <- subset_samples(ps, Organ=="Whole")
ps.filter.whole <- prune_taxa(taxa_sums(ps.filter.whole) >= 1, ps.filter.whole)
ps.filter.whole <- prune_samples(sample_sums(ps.filter.whole) >= 1, ps.filter.whole)
ps.filter.whole
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 7 taxa and 65 samples ]
## sample_data() Sample Data:       [ 65 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 7 taxa by 1 taxonomic ranks ]
# data pour plot
data_for_plot2 <- taxo_data_fast(ps.filter.whole, method = "abundance")
## Warning in `[<-.factor`(`*tmp*`, ri, value = "Other"): invalid factor level, NA
## generated
paste0("\n15 MOST ABUNDANT GENUS: \n") %>% cat()
## 
## 15 MOST ABUNDANT GENUS:
paste0("\"", levels(data_for_plot2$Name), "\",\n") %>% cat()
## "oligotype_AC (80) | size:2504 / N0253 (80) | size:2481.",
##  "oligotype_AG (109) | size:42030 / N0245 (109) | size:41133.",
##  "oligotype_AT (120) | size:3890567 / N0711 (120) | size:3744518.",
##  "oligotype_CC (92) | size:7065 / N0026 (92) | size:7069.",
##  "oligotype_GG (44) | size:4080 / N0250 (44) | size:4011.",
##  "oligotype_GT (111) | size:119651 / N0241 (111) | size:114576.",
##  "oligotype_TT (89) | size:15422 / N0244 (89) | size:15236.",
##  "Other.",
new_names <- c( "oligotype_AT (120) | size:3890567 / N0711 (120) | size:3744518.",
                "oligotype_GT (111) | size:119651 / N0241 (111) | size:114576.",
                "oligotype_AG (109) | size:42030 / N0245 (109) | size:41133.",
                "oligotype_TT (89) | size:15422 / N0244 (89) | size:15236.",
                "oligotype_CC (92) | size:7065 / N0026 (92) | size:7069.",
                "oligotype_GG (44) | size:4080 / N0250 (44) | size:4011.",
                "oligotype_AC (80) | size:2504 / N0253 (80) | size:2481."
)

data_for_plot2$Name <- factor(data_for_plot2$Name, levels = new_names)

col <- c("oligotype_AT (120) | size:3890567 / N0711 (120) | size:3744518."="#FEB24C", 
         "oligotype_GT (111) | size:119651 / N0241 (111) | size:114576."="#FAD769", 
         "oligotype_AG (109) | size:42030 / N0245 (109) | size:41133."="#666666", 
         "oligotype_CC (92) | size:7065 / N0026 (92) | size:7069."="#BEAED4", 
         "oligotype_TT (89) | size:15422 / N0244 (89) | size:15236."="#BF5B17",
         "oligotype_AC (80) | size:2504 / N0253 (80) | size:2481."="#F4CAE4", 
         "oligotype_GG (44) | size:4080 / N0250 (44) | size:4011."="#FDCDAC")

levels(data_for_plot2$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(data_for_plot2$Location) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))

data_for_plot2 <- data_for_plot2 %>% na.omit()

p2 <- ggplot(data_for_plot2, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, location=Location))+ 
  geom_bar(position = "stack", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text = element_text(size=14),
        #legend.key.height = unit(1, 'cm'),
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Location+Organ, scales = "free", ncol=3, labeller=label_parsed)+
  labs(x="Sample", y="Relative abundance", fill="Oligotype / MED node")

p2

Ovary (the most abundant nodes)

# select ovary
ps.filter.ovary <- subset_samples(ps, Organ=="Ovary")
ps.filter.ovary <- prune_taxa(taxa_sums(ps.filter.ovary) >= 1, ps.filter.ovary)
ps.filter.ovary <- prune_samples(sample_sums(ps.filter.ovary) >= 1, ps.filter.ovary)
ps.filter.ovary
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 7 taxa and 46 samples ]
## sample_data() Sample Data:       [ 46 samples by 18 sample variables ]
## tax_table()   Taxonomy Table:    [ 7 taxa by 1 taxonomic ranks ]
# data pour plot
data_for_plot3 <- taxo_data(ps.filter.ovary, method = "abundance", top=10, other=FALSE)
paste0("\n15 MOST ABUNDANT GENUS: \n") %>% cat()
## 
## 15 MOST ABUNDANT GENUS:
paste0("\"", levels(data_for_plot3$Name), "\",\n") %>% cat()
## "oligotype_AC (80) | size:2504 / N0253 (80) | size:2481.",
##  "oligotype_AG (109) | size:42030 / N0245 (109) | size:41133.",
##  "oligotype_AT (120) | size:3890567 / N0711 (120) | size:3744518.",
##  "oligotype_CC (92) | size:7065 / N0026 (92) | size:7069.",
##  "oligotype_GG (44) | size:4080 / N0250 (44) | size:4011.",
##  "oligotype_GT (111) | size:119651 / N0241 (111) | size:114576.",
##  "oligotype_TT (89) | size:15422 / N0244 (89) | size:15236.",
##  "Other",
new_names <- c( "oligotype_AT (120) | size:3890567 / N0711 (120) | size:3744518.",
                "oligotype_GT (111) | size:119651 / N0241 (111) | size:114576.",
                "oligotype_AG (109) | size:42030 / N0245 (109) | size:41133.",
                "oligotype_TT (89) | size:15422 / N0244 (89) | size:15236.",
                "oligotype_CC (92) | size:7065 / N0026 (92) | size:7069.",
                "oligotype_GG (44) | size:4080 / N0250 (44) | size:4011.",
                "oligotype_AC (80) | size:2504 / N0253 (80) | size:2481."
)

data_for_plot3$Name <- factor(data_for_plot3$Name, levels = new_names)

levels(data_for_plot3$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(data_for_plot3$Location) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))

data_for_plot3 <- data_for_plot3 %>% na.omit()

p3 <- ggplot(data_for_plot3, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, location=Location))+ 
  geom_bar(position = "stack", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text = element_text(size=14),
        #legend.key.height = unit(1, 'cm'),
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Location+Organ, scales = "free", ncol=3, labeller = label_parsed)+
  labs(x="Sample", y="Relative abundance", fill="Oligotype / MED node")

p3

Panels taxonomy of whole / ovary

legend_plot <- get_legend(p2 + theme(legend.position="bottom"))

# panels
p_group <- plot_grid(p2+theme(legend.position="none"), 
          p3+theme(legend.position="none"), 
          nrow=2, 
          ncol=1)+
    draw_plot_label(c("A1", "A2"), c(0, 0), c(1, 0.5), size = 20)

p_taxo <- plot_grid(p_group, legend_plot, nrow=2, rel_heights = c(1, .1))
p_taxo

Save taxonomic plot

setwd(path_plot)

tiff("2Da_OLIGO_counts_wolbachia.tiff", units="in", width=20, height=18, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
tiff("2Da_OLIGO_taxonomic_wolbachia_whole.tiff", units="in", width=16, height=12, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
tiff("2Da_OLIGO_taxonomic_wolbachia_ovary.tiff", units="in", width=18, height=14, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
tiff("2Da_OLIGO_taxonomic_wolbachia.tiff", units="in", width=18, height=16, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2
png("2Da_OLIGO_counts_wolbachia_big.png", units="in", width=20, height=18, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
png("2Da_OLIGO_counts_wolbachia_small.png", units="in", width=18, height=14, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
png("2Da_OLIGO_taxonomic_wolbachia_whole.png", units="in", width=16, height=12, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
png("2Da_OLIGO_taxonomic_wolbachia_ovary.png", units="in", width=18, height=14, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
png("2Da_OLIGO_taxonomic_wolbachia_big.png", units="in", width=18, height=18, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2
png("2Da_OLIGO_taxonomic_wolbachia_small.png", units="in", width=18, height=14, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2

Make main plot

setwd(paste0(path_oligo,"/wolbachia/oligotyping_Wolbachia_sequences-c2-s1-a0.0-A0-M10/HTML-OUTPUT"))

#setwd("/Users/hschrieke/Desktop/capture_ecran")
img <- magick::image_read("entropy.png")
p_entropy <- magick::image_ggplot(img, interpolate = TRUE)
p_entropy+ theme(plot.margin = unit(c(-7,-2.5,-7,-0.5), "cm"))

p_entropy+ theme(plot.margin=unit(c(-7,-2,-12,-5), "mm"))

aligned <- plot_grid(p_taxo, 
                     p_counts, 
                     align="hv")

aligned

p_entropy2 <- plot_grid(p_entropy, nrow=1)+
  draw_plot_label(c("C"), c(0), c(1), size=20, hjust=-0.5)

p_entropy2

t_plot <- plot_grid(aligned, 
                    p_entropy2,
                    nrow=2, 
                    ncol=1, 
                    scale=1,
                    rel_heights=c(2,1))

t_plot

setwd(path_plot)

tiff("2Da_OLIGO_main_wolbachia.tiff", width=36, height=36, res=300, units="in")
t_plot
dev.off()
## quartz_off_screen 
##                 2
png("2Da_OLIGO_main_wolbachia.png", width=36, height=36, res=300, units="in")
t_plot
dev.off()
## quartz_off_screen 
##                 2